Switching rates of multi-step reactions 
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We consider a switching rate of a meta-stable reaction scheme, which includes reactions with 
arbitrary steps, e.g. kA — > (k + r)A. Employing WKB approximation, controlled by a large system 
size, we evaluate both the exponent and the pre-exponential factor for the rate. The results are 
illustrated on a number of examples. 
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Since the celebrated work of Kramers on reaction-rate 
theory much effort has been devoted to extending 
and generalizing his results, see Ref. @ for a review. 
Applications of this theory can be found in fields as di- 
verse as high energy physics, nucleation, chemical kinet- 
ics, electric transport, diffusion in solids and population 
dynamics among many others. In this work we consider a 
switching rate in a generic reaction scheme, which admits 
more than one (quasi)stationary state. 

A particular case of single-step reactions allows for an 
exact solution and is well-studied in the literature [H, 0] • 
We thus concentrate on generic multi-step reactions. Al- 
though an exact solution is not known, a substantial 
progress may be achieved by adopting an analog of the 
quantum mechanical WKB approximation 0, 0, 0] , con- 
trolled by a large system size. With an exponential ac- 
curacy it gives the switching rate as an exponentiated 
action of an auxiliary mechanical problem. Evaluation 
of the pre-exponential factor requires a matching of the 
quasi-stationary distribution (QSD) function, found in 
the WKB framework, with the constant current "behind 
the barrier" solution [l], §]. The first consistent appli- 
cation of this strategy to a model reaction scheme was 
presented recently by Meerson and Sasorov Here we 
generalize their approach to an arbitrary scheme with 
metastable states. 

Consider a generic multi-step reaction scheme, where 
a state with n particles may be transformed into a state 
with n + r particles with the rate W r (n) . Here r is a set 
of integers not necessary equal ±1. The corresponding 
Master equation for the probability distribution P n (t) is 

d t P n (t) = J2[W r {n-r)P n - r (t)-W r (n)P n (t)} 

r 

= 52(e- rd "-l)W r (n)P n (t). (1) 

r 

We focus on reactions which admit a QSD centered at 
n = no and an unstable equilibrium (saddle point) at n = 
n s . For definiteness we assume that no < n s . We also 
assume that both no and n s scale in the same way with 
a large parameter N 3> 1, hereafter referred to as the 
system size, i.e. no, s ~ N. It is then convenient to pass 



to a scaling variable q — n/N and separate the leading 
and the first subleading orders in N in the corresponding 
reaction rates 

W r (n) = Nw r (q) + u r {q) + 0(l/N); q = n/N. (2) 

We seek for QSD in the form P n (t) = ir(n)e ot , where 
Eo = 1/r is an exponentially small eigenvalue of the Mas- 
ter equation. In the rescaled coordinate the correspond- 
ing eigenvector may be sought in the WKB form 



n(q) =cxp{-NS(q)-S 1 (q)} 



(3) 



Substituting this form in the Master equation ([T]) and 
keeping terms up to the order of 0(1), one finds 



rS' 



U r ) 



2N b + N bl N w r 



1 , (4) 



where the primes denote derivatives with respect to 
rescaled reaction coordinate q. We have also took into 
account that the eigenvalue E is expected to be expo- 
nentially small in N (see below) and thus may be omitted. 

In the order N this equation acquires a form of the 
stationary Hamilton-Jacobi equation H (q, S') — 0, where 
the effective classical Hamiltonian takes the form [a, M 



fl-(g,p) = X>r(g) (e rp -i) 



(5) 



and we have denoted S' — p. Therefore to the order N 
the problem is reduced to finding zero energy trajectories 
p = p(q), such that H(q,p(q)) = 0, of a corresponding 
" mechanical" problem. 

The phase portrait of a typical bistable reaction is plot- 
ted in Fig. [TJ There are at least two appropriate zero 
energy trajectories: the relaxation trajectory p = and 
the activation trajectory p = p a (q), see Fig. [TJ The 
classical equation of motion along the relaxation path 
q = H p (q,0) = ^2 r rw r (q) is nothing but the mean-field 
rate equation for our reaction scheme. According to our 
assumptions it admits stationary states qo, s — no,s/N, 
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FIG. 1: (Color online) Phase portrait of a typical bistable 
reaction. The dashed line is the relaxation trajectory p = 0, 
the solid line is the activation trajectory p a = p a (q)- These 
zero-energy curves intersect at the metastable points qo, q' 
and the saddle point q 3 . The arrows show direction of motion 
according to the Hamilton equations. The matching of acti- 
vation and relaxation solutions takes place in a narrow region 
of the width l a ~ iV" 1 / 2 around the saddle point. 



where H p (qo^ s ,0) = (other stationary states are pos- 
sible, e.g. q' see Fig. Q]). Those are the points, where 
the activation trajectory p a (q) crosses the relaxation one 
p = and thus p a {lo,s) = 0. 

To escape from a metastable state centered around qq 
the system must evolve along the activation trajectory, 
Fig. Q] The QSD is given by Eq. ©, where S(q) and 
Si(q) are determined by the order TV and order 1 terms 
in Eq. Q correspondingly. They lead to 



S(q) 
Si(q) 



dqp a {q) 



dq '■ 



p' a H pp + 2H pq -2Y, r u r {e r 



1) 



2H„ 



(0) 
(7) 



where derivatives of the Hamiltonian are evaluated along 
the activation path, e.g. H pq — ^2 r re rPa ^w' r (q), etc 
and p' a = S". Equations ©, ©, determine QSD up 
to a multiplicative constant. To find the latter, one needs 
to match the QSD with the constant current solution on 
the other side of the saddle point q s 

At q > q s the system evolves along the relaxation tra- 
jectory p — 0, Fig. [I] and therefore S = 0. Solving Eq. 
for Si, one finds 



n(q) = J/H p (q,0), 



(8) 



where J is an integration constant given by the current 
out of QSD. Indeed, the Master equation (fTJ) , having the 
structure of the continuity relation, in a vicinity of the 
relaxation trajectory p = acquires a form 

d t P(q, t) = -d q [H p (q, 0)P(q, t) + 0(l/N)] . (9) 



Therefore the relaxation limit © of QSD P(q,t) = 
Tr(q)e~ Eot is nothing but a constant current, J, solution 
of the Master equation (where we have again neglected 
the exponentially small eigenvalue Eq on the l.h.s.). On 
the other hand, integrating the continuity relation ^ 
throughout the region of support of QSD and assuming 
that escape takes place only through the saddle point q s 
[l(|, one finds 



E J-K(q)dq = J . 



(10) 



Finally to establish relation between the activation so- 
lution, Eqs. (|3|), (O, ([7|), at q < q s and the relaxation 
one, Eq. (|8|), at q > q s , one needs to consider Master 
equation in an immediate vicinity of the saddle q s 0]. 
Expanding the r.h.s. of Eq. (JTJ) to the second derivative, 
one finds for the current: 

[H pq (q s , 0)] (q - q s )n(q) - Hp jp ° } d g n(q) = J , (11) 

where we have used the fact that at the saddle point 
H p (q s ,0) = J]r™r(<7s) = °- Solution of Eq. ([II]) with 
a proper asymptotic behavior has the following form 
tt(?) = {2NJ/H pp )e^-"^ 2 / 1 ' f™ qs dqe-to-^'ft, where 
if = H pp (q s ,0)/NH pq (q s ,0). Indeed, away from the sad- 
dle point q s it possesses the following asymptotics: 



ir{q) = 



{q-q,)H pq ' 



2NJl 3 V^ p (q-q s ) 2 /l 2 s 



q- <? s > l s , 



q^>h 



(12) 



The first line matches with the relaxation solution (|8j at 
q w q s , as it should. The second line is to be matched 
with the activation solution Eqs. ([3]), ©, which in 
the vicinity of q = q s takes the form 



n(q) = e - NSl - q ^- Sl( - g ^ e -N{q-q s fS"(q B )/2^ 



(13) 



To relate the g-dependent exponential factors here and 
in the second line of Eq. fT2"]) one may differentiate the 
identity H(q,p a (q)) = over q to find 



H, 



!~rHpp' a 



0: 



H qq + E p p" a + (H pp p a + 2H pq )p' a = . 

(14) 

Employing that p' a = S" and H(q,Q) = H p (q o ^ s ,0) = 0, 
one finds 



2H pq (q ^ s ,0) 2Y,r rw 'r{qo,s) 
* (Qo,s) — — - 



Hpp(qo, s ,0) 



J2 r r 2 w r (q , s ) 



(15) 



and therefore S"(q s ) = —2/Nl 2 s . This establishes equal- 
ity of the exponential factors in Eqs. (fT2)) and (fT3|) . Com- 
paring the pre-exponential coefficients one finds for the 
escape current: 



Hpp(qs,0) l \S^qs)\ e -iVS fe )- Sl ( 9s ) _ (16) 
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One can employ now the normalization condition ([10]) 
to find the escape rate Eq — 1/r. To this end we notice 
that the bulk of the QSD is centered around q and ap- 
proximate the integral by the Gaussian one. As a result 
one finds for the escape time 



4tt 



oSl(q s )-S 1 (q ) 



,N[S(q s )-S(q )] 



H pP (q s ,0) y/\S"(q a )\S"(q ) 



, (17) 



where S(q s ) — S(qo) and Si(q s ) — Si(qo) are fully deter- 
mined by Eqs. (J6j> and ([7|). It is important to mention 
that the corresponding integrals are free of singularities 
and can be straightforwardly evaluated for any given re- 
action scheme. Equation (fT"7|) is a main result of the 
present letter. 

For analytically tractable examples it is useful to notice 
that, with the help of identities (jT4j) one may partially 
integrate Eq. ([7]) to obtain 



Si(q s ) - S 1 (q ) 



In . 



S"(qo) 
\S"(q s )\ 



(18) 



dq 



H qq 
2H„ 



1) 



Employing Eq. (fT5|) , one may somewhat simplify Eq. (fT7|) 
to cast it in the following form 



2?re 



A 



H pq (q s ,0) 



D N[S(q s )-S(q )] 



(19) 



Below we illustrate usefulness of Eqs. (fT"7|) and (fl"9|) on a 
few examples. 

r\^T2 reactions. Consider a reaction scheme, where 
the step variable r may acquire only two values r\ and 
r2 . The corresponding reaction rates are W ri 2 (n) = 
Nw ri2 (q), where we have omitted possible subleading 
terms u ri 2 for brevity. The Hamiltonian takes the form 

H(q, p) = w ri (q) (e™ -l) + w r2 (?) (e™ - 1) , (20) 

and the activation trajectory is given by the solution of 
the following algebraic equation for e Pa 



e riPa{q) _ I 

6 r2Pa(q) — 1 



Wr 2 (q) 
w ri [q) ' 



(21) 



As a result, the following identity holds along the activa- 
tion trajectory: 



<i<i 



<(g)(e rip ° -l) + <(g)(e^° -1) 
w'(q)(e^P« - l) + w>(q)(e r *P° - 1) 



= — ln{w ri w - w w r2 ) 

w ri w' r2 — w' ri w r2 dq 



The fixed points satisfy: riw ri (qo iS ) — — T2W T2 (qo,s), 
while H pq (q QtS ,0) = rxw'^iqo^) + r 2 w' r2 (q 0iS ). Employing 
Eqs. (fT8| and (fT9|) . one finds for the switching time 



T = 2lT\ 



(Qo) 



,N[S(q s )-S{q )] 



^\H pq (q s ,0)H pq (q o ,0)\ ' 



(22) 



where w ri (q s )/w ri (qo) = w r2 (q s )/w r2 (qo) and the action 
is given by Eq. (J6]). 

In a particular case of single-step reactions, = ±1, 
Eq. (f2"Tj) may be solved explicitly, e Pa ^ = W-(q)/w + (q). 
The fixed points are given by w+(qo. s ) = tu_(qo jS ) and 
according to Eq. ([T5]) H pq (q o , s ,0) = -S"(qo, s )w + (qo, s )- 
Employing Eq. (|22|) . the switching rate of the single-step 
reaction schemes may be written as 



-/£*r(^-^) 



W±(q ) y/\S"(q a )\S"(q ) 



2tt 



o N[S(q B )-S{q )] 



(23) 



where 



S(q s )-S(q )= dqln(w-(q)/w + (q)) (24) 
J qo 

and we have included subleading terms in the rates u± (q) , 
according to Eq. (fT8|) , [ll[ . In a particular case of reac- 
tion rates having only leading terms (u± = 0) Eq. |23|) 
coincides with results of Doering et al. [4], who have 
shown it to be the large N asymptotic of the exact result 
for the single-step reactions [3j. In general, the u r terms 
can substantially modify the prefactor [§] (see below). 
Demographic explosion. Consider a single-step model 
A <=» with the relative rates 1 and AT(1 - <5 2 )/2, 
where < S < 1, and 2 A — > 3 A with the relative rate 
1/N. The corresponding transition rates are 



W-{n)=n- W+{n) = 



N(l-5 2 ) n(n-l) 
2 + 2N 



The rescaled rates are W- = q; w + = (1 — 5 2 + q 2 )/2 
while it_ — and u + = —q/2 and the two rescaled fixed 
points are qo, s — 1 T 8. Employing Eq. (|23|) . one finds 
for the escape time from the metastable state centered 
at n = N(l — S) towards n — > oo 



2tt 1 



= JV[S(l+(5)-S(l-*)] 



S 1 



(25) 



in a perfect agreement with Meerson and Sasorov recent 
result 0] . This example is specially interesting because it 
shows the importance of the subleading terms u r . Disre- 
garding these terms, one obtains a prefactor proportional 
to (1 — (J) -1 / 2 instead of the correct one (1 — S)^ 1 . This 
constitutes an arbitrarily large error in the limit 5 — > 1, 
where the action 5(2) — 5(0) remains bounded. 

Fokker-Planck Hamiltonian. Consider a dissipative 
particle under an influence of a multiplicative Gaussian 
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noise (understood in the sense of Ito Q). The corre- 
sponding Fokker-Planck equation is P = HP, where 



H(q,p) =p 2 D(q)-pV'(q), 



(26) 



here D(q) > is a coordinate-dependent diffusion coef- 
ficient and p = —d q . Since this is a normally ordered 
operator, cf. Eq. ([T]), one may employ the theory devel- 
oped above. Following WKB approximation one substi- 
tutes p — > p and employs Eq. (| 1 9|1 . The stationary points 
are defined by the condition V'(qo, s ) = and the acti- 
vation trajectory is given by p a (q) — V'(q)/D(q). As a 
result S(q s )-S(q ) = QdqV'(q)/D(q) and H pq (q s , 0) = 
—V"(q s ) > 0. There are no subleading terms here, 
u r = 0, and therefore 



A = 



qs H 1 



V"(q s )D(q s ) 



V"(q )D( qo ) 



where we have made use of V'(qo) — V'(q s ) = 0. Using 
this equality again one finds S"(qo, s ) — U"(go,s)/-D(<7o,s), 
and finally, plugging all together in Eq. f(T9j) , one obtains 



2tt 



! D(q s 



,/V"(qo)\V"{q s )\ V D (lo) 



e Jl'dqV'(q)/D(q) 



(27) 



in agreement with previous calculations Q. Assuming 
a constant diffusion coefficient D(q) = T (i.e. additive 
noise), one recovers Kramers result Notice that the 
role of N is played by 1/T. 

Higher moments of noise. Consider now Kramers 
problem of a dissipative particle subject to a white, non- 
Gaussian noise. The corresponding Hamiltonian reads 
as 



H(q,p)=e k p k + Tp 2 -pV'(q) 



(28) 



Here k = 3,4, • • • and £3,4,... is the third, fourth, etc (i.e. 
first non-vanishing beyond the second) irreducible mo- 
ment of the noise correlation function. This type of noise 
appears as e.g. higher order corrections in the Kramers- 
Moyal expansion of the master equation Assuming 
that the higher moments are much smaller than the sec- 
ond one [l2j and proceeding as in the last case we find 



2tt 



y/\V"(q s )\V"(q Q ) 



exp 



,(V(qs)-V(qo))/T 



[V\q)f- 1 dq + 0(et) 



'in 



(29) 



As can be seen, the prefactor remains unchanged and the 
whole contribution coming from the higher order noise 
concentrates in an extra " phase" . Note that is neces- 
sarily positive for even k (in order to keep the noise real) 
but it can be either positive or negative for odd k. For 



the escape processes under consideration V(q s ) > V(qo), 
and so the integral term in the extra "phase" is positive, 
what implies that even moments of noises only contribute 
to reduce the escape time, while the odd ones can reduce 
or increase the switching time, depending on the sign of 
the corresponding moment. 

To conclude we have calculated the escape rate from 
a metastable state whose dynamics is described by a 
general multi-step master equation. We found a rela- 
tively simple analytical result for switching rates between 
metastable states (but not for absorbing phase transition, 
as e.g. extinction) of an arbitrary single-species reac- 
tion scheme. We have shown that the general formula 
found here reduces to known results for single-step reac- 
tions and Fokker-Planck equations, with either additive 
or multiplicative noises. 

We are indebted to M. Dykman, B. Meerson and P. 
Sasorov for numerous useful discussions. C. E. is grateful 
to the William I. Fine Theoretical Physics Institute for its 
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A.K. was supported by NSF grants DMR-0405212 and 
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In so doing we neglect the back-current from the state 
centered at q' , see Fig.[T] towards qo. This is appropriate 
at times shorter than (exponentialy long) escape time 
from qo- 

Notice that equation of motion along the activation tra- 
jectory q = Hp(q,p a (q)) = W- — w + = —H p (q,0) is the 
time-reversal partner of the relaxation motion. This is 
a consequence of the fact that the single-step reactions 
satisfy detailed balance condition. We are indebted to M. 
Dykman for discussion of this point. 
Here, weakness is to be understood in the sense that 
0(e|) terms in the exponent may be disregarded. This 
is the case when t\{k - 



\)f qs [V'(q)] 2k ~ 3 dq<^T 2k -' 



